from Algo_var import *
import work_var
from TimingLib import TimingLib
import numpy as np
import math
import pdb
import warnings
from ctarnoldi import build_matrix_dr, ctarnoldi, compute_poles_res

"""
modified from https://github.com/gzz2000/DelayModelingCollections
simple LUT class with linear interpolation
"""
class LUT:
    def __init__(self, xs, ys, vs):
        self.xs = xs
        self.ys = ys
        self.vs = vs
        assert len(self.vs) == len(self.xs) * len(self.ys)

    def transpose(self):
        xs, ys = self.ys, self.xs
        vs = [0. for i in range(len(self.vs))]
        for i in range(len(xs)):
            for j in range(len(ys)):
                vs[i * len(ys) + j] = self.vs[j * len(xs) + i]
        return LUT(xs=xs, ys=ys, vs=vs)

    def lookup_grad(self, x, y):
        n = len(self.xs)
        m = len(self.ys)
        xi, xj, yi, yj = 0, 0, 0, 0
        if n > 1:
            while xi + 1 < n and self.xs[xi + 1] <= x: xi += 1
            xj = xi
            if xj + 1 < n: xj += 1
            else: xi -= 1
        if m > 1:
            while yi + 1 < m and self.ys[yi + 1] <= y: yi += 1
            yj = yi
            if yj + 1 < n: yj += 1
            else: yi -= 1
        def FT(i, j):
            return self.vs[i * m + j]
        def XT(i):
            return self.xs[i]
        def YT(j):
            return self.ys[j]
        rxi = interpolate(FT(xi, yi), FT(xi, yj), YT(yi), YT(yj), y)
        rxj = interpolate(FT(xj, yi), FT(xj, yj), YT(yi), YT(yj), y)
        ryi = interpolate(FT(xi, yi), FT(xj, yi), XT(xi), XT(xj), x)
        ryj = interpolate(FT(xi, yj), FT(xj, yj), XT(xi), XT(xj), x)
        r = interpolate(rxi, rxj, XT(xi), XT(xj), x)
        if abs(rxi - rxj) < 1e-10:
            dx = 0
        else:
            dx = (rxj - rxi) / (XT(xj) - XT(xi))
        if abs(ryi - ryj) < 1e-10:
            dy = 0
        else:
            dy = (ryj - ryi) / (YT(yj) - YT(yi))
        return r, dx, dy

    def lookup(self, x, y):
        return self.lookup_grad(x, y)[0]

def interpolate(f1, f2, x1, x2, x):
    if x1 == x2: return f1
    else: return ((x2 - x) * f1 + (x - x1) * f2) / (x2 - x1)

def apply_derate(vs, derate):
    return [v * derate for v in vs]

class Delay_Calculator:
    """Ceff, delay and transition calculation"""
    """
    self.design -> design name of current design
    self.SpefNets -> spef nets of current design
    self.Parasitics -> rctree parastic class for delay calculation
    self.rc_corner -> RC corner of spef
    self.th_slew1 -> lower threshold of slew
    self.th_slew2 -> upper threshold of slew
    self.th_delay -> threshold of delay
    self.Elmore -> calculated Elmore delay
    
    function:
    self.LUT_generate: generate LUT class based on Delay or Trans
    self.Calcu_Pin_Ceff_Delay_Trans: core algorith using Newton regression for ceff, delay, slew calculation
    self.calc_waveform: calculate waveform of certain node on RC tree
    self.calc_waveform_grad: calculate grad of waveform
    self.calc_fval: residual function for Newton regression calculation (J(xk) * delat xk = -f(xk))
    self.calc_jacobian: calculate jacobian matrix J(x)
    self.test_Ceff_smallNR and test_Ceff_ultrasmallNR: different NR with 3 or 2 equations
    self.calc_single_timepoint_nr: calculate timing point of wave at certain voltage using NR
    self.Calc_Elmore: calculate Elmore delay of a given net
    """
    def __init__(
        self,
        design: str = "mc_top",
        rc_corner : str = "rcworst_CCworst",
        th_slew1 : float = 0.3,
        th_slew2: float = 0.7,
        th_delay: float = 0.5,
        ):
        print(f"Reading spef nets of design {design}")
        self.design = design
        self.rc_corner = rc_corner
        Spef = os.path.join(work_var.inn_data_dir, f"{design}/{design}_postRoute_fill_std_{rc_corner}.spef")
        self.Parasitics = Spef_Parser.Read_Spef(Spef)
        self.th_slew1 = th_slew1
        self.th_slew2 = th_slew2
        self.th_delay = th_delay
        self.Elmore = {}
    
    def LUT_generate(self, NLDM):
        lut = [] #cell delay or trans   
        if isinstance(NLDM, TimingLib_Parser.Delay):
            for i in NLDM.delay:
                for j in i:
                    lut.append(j)
        elif isinstance(NLDM, TimingLib_Parser.Trans):
            if NLDM.isscalar:
                raise ValueError("Trans object with is_scalar=True is not allowed.")
            else:
                for i in NLDM.trans:
                    for j in i:
                        lut.append(j)
        else:
            raise TypeError(f"Expected Delay or Trans, but got {type(NLDM)}")
        
        return LUT(NLDM.index1, NLDM.index2, lut).transpose()    
    
    # Gate delay computation with CT-Arnoldi reduced RC loads
    # adapted from Dartu, Menezes, Pileggi, Performance Computation for Prec aracterized CMOS Gates with RC Loads, TCAD 1996 
    
    def Calcu_Pin_Ceff_Delay_Trans(self, parasitic, input_slew, Delay_LUT, Trans_LUT, order=1, cal_net=True):
        tree_rc = parasitic
        
        # initialize parameters
        total_cap = sum(cap for _, cap in tree_rc.grounded_caps) + sum(tree_rc.sink_cell_caps)
        delay_lut = self.LUT_generate(Delay_LUT)
        slew_lut = self.LUT_generate(Trans_LUT)
        
        # adapted from eq(22) to use slew instead.
        """
        Fitting voltage model v(t) = V0 * e^(-t / Rd * CL)
        -> ln(v(t)/V0) = -t/Rd*CL -> ln(th1) = -t1/Rd*CL, ln(th2) = -t2/Rd*CL
        -> slew = t2 - t1 = -Rd * CL * ln(th2/th1)
        -> 0.9 is a empirical coefficient
        """
        init_cap = min(total_cap, delay_lut.xs[-1])
        init_slew = slew_lut.lookup(init_cap, input_slew)
        init_delay = delay_lut.lookup(init_cap, input_slew)
        Rd = init_slew / (init_cap * math.log(self.th_slew2 / self.th_slew1)) * 0.9
        C, G = build_matrix_dr(tree_rc, Rd)
        # reduction of C, G
        Uq, Hq, Glu, Gpiv = ctarnoldi(C, G, order, Rd)
        # calculate poles and residues
        flag_value = False
        #ugly fix to value error
        try:
            poles, residues_mat = compute_poles_res(Uq, Hq, C, G, Glu, Gpiv, Rd)
        except ValueError as e:
            print(f"Captured exception: {e}, type={type(e)}")
            flag_value = True
        # adapted from eq(16)-(17). initial guess of Ceff and driver
        Ceff = init_cap
        dt = init_slew / (self.th_slew2 - self.th_slew1)
        # t0: Gate intrinsic delay
        t0 = init_delay - 0.69 * Rd * Ceff - dt / 2.

        # our main task is to solve the following equation of
        #  Ceff, t0, and dt.  (adapted from eq(12))
        # Ipi(dt, dt) = ICeff(dt, dt, Ceff) ..... [0]
        # Timeof(0.5, t0, dt) = delay(Ceff)
        #  -> i.e., waveform(t0, dt, delay(Ceff)) = 0.5 ..... [1]
        # Timeof(0.8, t0, dt) - Timeof(0.2, t0, dt) = slew(Ceff)
        #
        # introduce an auxiliary variable tr1,
        # and replace the last equality with the following:
        # waveform(t0, dt, tr1) = 0.2 ..... [2]
        # waveform(t0, dt, tr1 + slew(Ceff)) = 0.8 ..... [3]
        tr1 = init_delay - init_slew * (1 - self.th_slew1)
        if not flag_value:
            # ugly fix to ensure no overflowError occur
            flag_overflow = False
            try:
                Ceff, delay, slew, dt, t0 = self.test_Ceff_ultrasmallNR(Ceff, dt, t0, Rd, residues_mat, poles, input_slew, delay_lut, slew_lut)
            # ceff is too small when overflow occur
            except OverflowError:
                Ceff = total_cap
                delay = delay_lut.lookup(Ceff, input_slew)
                slew = slew_lut.lookup(Ceff, input_slew)
                flag_overflow = True
            #print(Ceff,delay,slew)
            
            if cal_net:
                fanout_slew, fanout_delay = {}, {}
                for endpoint, io in tree_rc.endpoints:
                    if io != 'I': continue
                    if not flag_overflow:
                        try:
                            t_delay = self.calc_single_timepoint_nr(dt, endpoint, self.th_delay, delay - t0, residues_mat, poles)
                            t_slew1 = self.calc_single_timepoint_nr(dt, endpoint, self.th_slew1, delay - t0, residues_mat, poles)
                            t_slew2 = self.calc_single_timepoint_nr(dt, endpoint, self.th_slew2, delay - t0, residues_mat, poles)
                            # output slew of endpoint
                            fanout_slew[tree_rc.names[endpoint]] = t_slew2 - t_slew1
                            # net arc delay, 0.7 is a coefficient to mitigate net delay error
                            fanout_delay[tree_rc.names[endpoint]] = (t_delay + t0 - delay) * 0.7
                            #print(f'Node {tree_rc.names[endpoint]} ({endpoint}): delay={t_delay + t0}, slew={t_slew2 - t_slew1}')
                        except ZeroDivisionError:
                            fanout_slew[tree_rc.names[endpoint]] = slew
                            fanout_delay[tree_rc.names[endpoint]] = delay
                    else:
                        fanout_slew[tree_rc.names[endpoint]] = slew
                        fanout_delay[tree_rc.names[endpoint]] = delay
            else:
                fanout_slew, fanout_delay = None, None
        else:
            Ceff = total_cap
            delay = delay_lut.lookup(Ceff, input_slew)
            slew = slew_lut.lookup(Ceff, input_slew)
            if cal_net:
                fanout_slew, fanout_delay = {}, {}
                for endpoint, io in tree_rc.endpoints:
                    if io != 'I':
                        continue
                    fanout_slew[tree_rc.names[endpoint]] = delay
                    fanout_delay[tree_rc.names[endpoint]] = slew
            else:
                fanout_slew, fanout_delay = None, None
        return Ceff, delay , slew , fanout_slew, fanout_delay
    
    def calc_waveform(self, t, dt, residues_mat, poles, node_id=0):
        if t < 0: return 0.
        t1 = max(0, t - dt)
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=RuntimeWarning)
            return (t - t1 - np.dot(residues_mat[node_id] / poles, np.exp(t * poles) - np.exp(t1 * poles))) / dt

    def calc_waveform_grad(self, t, dt, residues_mat, poles, node_id=0):
        if t < 0: return 0., 0.
        # some repeated calculation can be avoided here.
        g_dt = -self.calc_waveform(t, dt, residues_mat, poles, node_id=node_id) / dt
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=RuntimeWarning)
            g_t = (1. - np.dot(residues_mat[node_id], np.exp(np.clip(t * poles, -700, 700)))) / dt
        if t >= dt:
            g_left = (1. - np.dot(residues_mat[node_id], np.exp(np.clip(((t - dt) * poles), -700, 700)))) / dt
            g_t -= g_left
            g_dt += g_left
        return g_t, g_dt

    def calc_fval(self, Ceff, dt, t0, tr1, residues_mat, poles):
        fval = np.zeros(4, dtype=np.float32)
        
        Qpidt = np.dot(residues_mat[0], (-1. + np.exp(dt * poles) - dt * poles) / (dt * poles**2 * Rd))
        QCeffdt = Ceff + Rd * Ceff**2 * (-1. + math.exp(-dt / (Ceff * Rd))) / dt
        fval[0] = Qpidt - QCeffdt

        delay = delay_lut.lookup(Ceff, input_slew)
        slew = slew_lut.lookup(Ceff, input_slew)
        
        fval[1] = self.calc_waveform(delay - t0, dt, residues_mat, poles) - self.th_delay
        fval[2] = self.calc_waveform(tr1 - t0, dt, residues_mat, poles) - self.th_slew1
        fval[3] = self.calc_waveform(tr1 + slew - t0, dt, residues_mat, poles) - self.th_slew2

        return fval, delay, slew

    def calc_jacobian(self, Ceff, dt, t0, tr1, residues_mat, poles):
        delay, delay_grad_x, _ = delay_lut.lookup_grad(Ceff, input_slew)
        slew, slew_grad_x, _ = slew_lut.lookup_grad(Ceff, input_slew)

        jacobian = np.zeros((4, 4), dtype=np.float32)
        
        # df0 / dCeff = - dQCeffdt / dCeff
        jacobian[0, 0] = -(dt - 2. * Ceff * Rd + math.exp(-dt / (Ceff * Rd)) * (dt + 2 * Ceff * Rd) / dt)
        # df0 / ddt
        jacobian[0, 1] = (
            # dQpidt / ddt
            np.dot((1. + np.exp(dt * poles) * (-1. + dt * poles)) / (dt**2 * poles**2 * Rd), residues_mat[0]) -
            # dQCeffdt / ddt
            Ceff * (np.exp(-dt / (Ceff * Rd)) * (-dt - Ceff * Rd) + Ceff * Rd) / dt**2
        )

        # df{1, 2, 3}
        gwf1t, gwf1dt = self.calc_waveform_grad(delay - t0, dt, residues_mat, poles)
        gwf2t, gwf2dt = self.calc_waveform_grad(tr1 - t0, dt, residues_mat, poles)
        gwf3t, gwf3dt = self.calc_waveform_grad(tr1 + slew - t0, dt, residues_mat, poles)

        # df / ddt
        jacobian[1, 1] = gwf1dt
        jacobian[2, 1] = gwf2dt
        jacobian[3, 1] = gwf3dt

        # df / dt0
        jacobian[1, 2] = -gwf1t
        jacobian[2, 2] = -gwf2t
        jacobian[3, 2] = -gwf3t

        # df / dCeff
        jacobian[1, 0] = gwf1t * delay_grad_x
        jacobian[3, 0] = gwf3t * slew_grad_x

        # df / dtr1
        jacobian[2, 3] = gwf2t
        jacobian[3, 3] = gwf3t

        return jacobian

    def calc_f123val(self, delay, slew, dt, t0, tr1, residues_mat, poles):
        fval = np.zeros(3, dtype=np.float32)
        fval[0] = self.calc_waveform(delay - t0, dt, residues_mat, poles) - self.th_delay
        fval[1] = self.calc_waveform(tr1 - t0, dt, residues_mat, poles) - self.th_slew1
        fval[2] = self.calc_waveform(tr1 + slew - t0, dt, residues_mat, poles) - self.th_slew2
        return fval

    def calc_f123jacobian(self, delay, slew, dt, t0, tr1, residues_mat, poles):
        jacobian = np.zeros((3, 3), dtype=np.float32)
        
        # df{1, 2, 3}
        gwf1t, gwf1dt = self.calc_waveform_grad(delay - t0, dt, residues_mat, poles)
        gwf2t, gwf2dt = self.calc_waveform_grad(tr1 - t0, dt, residues_mat, poles)
        gwf3t, gwf3dt = self.calc_waveform_grad(tr1 + slew - t0, dt, residues_mat, poles)

        # df / ddt
        jacobian[0, 0] = gwf1dt
        jacobian[1, 0] = gwf2dt
        jacobian[2, 0] = gwf3dt

        # df / dt0
        jacobian[0, 1] = -gwf1t
        jacobian[1, 1] = -gwf2t
        jacobian[2, 1] = -gwf3t

        # df / dtr1
        jacobian[1, 2] = gwf2t
        jacobian[2, 2] = gwf3t

        return jacobian

    def calc_f12val(self, delay, slew, dt, t0, residues_mat, poles):
        fval = np.zeros(2, dtype=np.float32)
        fval[0] = self.calc_waveform(delay - t0, dt, residues_mat, poles) - self.th_delay
        fval[1] = self.calc_waveform(delay - slew * (self.th_delay - self.th_slew1) / (self.th_slew2 - self.th_slew1) - t0, dt, residues_mat, poles) - self.th_slew1
        return fval

    def calc_f12jacobian(self, delay, slew, dt, t0, residues_mat, poles):
        jacobian = np.zeros((2, 2), dtype=np.float32)
        gwf1t, gwf1dt = self.calc_waveform_grad(delay - t0, dt, residues_mat, poles)
        gwf2t, gwf2dt = self.calc_waveform_grad(delay - slew * (self.th_delay - self.th_slew1) / (self.th_slew2 - self.th_slew1) - t0, dt, residues_mat, poles)
        # df1 / d{dt, t0}
        jacobian[0, 0] = gwf1dt
        jacobian[0, 1] = -gwf1t
        # df2 / d{dt, t0}
        jacobian[1, 0] = gwf2dt
        jacobian[1, 1] = -gwf2t
        return jacobian

    def test_Ceff_smallNR(self, Ceff, dt, t0, tr1, input_slew, residues_mat, poles):
        delay_old = None  # detect convergence
        slew_old = None   # detect convergence
        Ceff_old = None
        dt_old, t0_old, tr1_old = None, None, None
        for ceff_i in range(200):
            delay = delay_lut.lookup(Ceff, input_slew)
            slew = slew_lut.lookup(Ceff, input_slew)
            # print('===============================')
            # print(f'Ceff iter {ceff_i}: Ceff={Ceff}, delay={delay}, slew={slew}')
            if delay_old is not None and \
            math.fabs((delay_old - delay) / delay_old) < 1e-3 and \
            math.fabs((slew_old - slew) / slew_old) < 1e-3:
                # print(f'| CONVERGED (1‰). STOP.')
                break

            nr_converged = False
            for nr_i in range(20):
                fval = self.calc_f123val(delay, slew, dt, t0, tr1, residues_mat, poles, delay_lut, slew_lut)
                # print(f'\\__ smallNR iter {nr_i}, dt={dt}, t0={t0}, tr1={tr1}, fval={fval}')
                if np.max(np.abs(fval)) < 1e-3:
                    # print(f'\\__ smallNR EARLY STOP.')
                    nr_converged = True
                    break
                jacobian = self.calc_f123jacobian(delay, slew, dt, t0, tr1, residues_mat, poles)
                delta = np.linalg.pinv(jacobian) @ fval
                dt -= delta[0]
                t0 -= delta[1]
                tr1 -= delta[2]
                # print(f'\\__     UPDATE delta={delta}')
                if dt < 0.:
                    # print('WARN: Too small dt, Rd too large or Ceff decreased too fast')
                    break
            if not nr_converged:
                # print('WARN: smallNR not converged after many iters.')
                if Ceff_old is not None:
                    # print('WARN: -> trying to lower Ceff step..')
                    Ceff = (Ceff + Ceff_old) / 2.
                    dt, t0, tr1 = dt_old, t0_old, tr1_old
                    continue
                
            Qpidt = np.dot(residues_mat[0], (-1. + np.exp(dt * poles) - dt * poles) / (dt * poles**2 * Rd))
            QCeff_quad_coeff = Rd * (-1. + math.exp(-dt / (Ceff * Rd))) / dt
            Ceff_new = (-1. + math.sqrt(1 + 4. * QCeff_quad_coeff * Qpidt)) / (2. * QCeff_quad_coeff)
            # print(f'| NEW Ceff={Ceff_new} (old: {Ceff})')
            # assert Ceff_new < Ceff, 'Should most likely be decreasing.'
            Ceff_old = Ceff
            Ceff = Ceff_new
            delay_old = delay
            slew_old = slew
            dt_old, t0_old, tr1_old = dt, t0, tr1
        else:
            print('WARN: Ceff and delay not converged after many iters.')

        return Ceff, delay, slew, dt, t0, tr1

    def test_Ceff_ultrasmallNR(self, Ceff, dt, t0, Rd, residues_mat, poles, input_slew, delay_lut, slew_lut):
        delay_old = None  # detect convergence
        slew_old = None   # detect convergence
        Ceff_old = None
        dt_old, t0_old = None, None

        # calculate the Q in [0, dt / REGION] instead of
        # [0, dt]. an ugly fix to the wrong Ceff..
        REGION = 0.55
        
        for ceff_i in range(200):
            delay = delay_lut.lookup(Ceff, input_slew)
            slew = slew_lut.lookup(Ceff, input_slew)
            # print('===============================')
            # print(f'Ceff iter {ceff_i}: Ceff={Ceff}, delay={delay}, slew={slew}')
            if delay_old is not None and \
            math.fabs((delay_old - delay) / delay_old) < 1e-3 and \
            math.fabs((slew_old - slew) / slew_old) < 1e-3:
                # print(f'| CONVERGED (1‰). STOP.')
                break

            nr_converged = False
            for nr_i in range(20):
                fval = self.calc_f12val(delay, slew, dt, t0, residues_mat, poles)
                # print(f'\\__ smallNR iter {nr_i}, dt={dt}, t0={t0}, fval={fval}')
                if np.max(np.abs(fval)) < 1e-3:
                    # print(f'\\__ smallNR EARLY STOP.')
                    nr_converged = True
                    break
                jacobian = self.calc_f12jacobian(delay, slew, dt, t0, residues_mat, poles)
                delta = np.linalg.pinv(jacobian) @ fval
                dt -= delta[0]
                t0 -= delta[1]
                # print(f'\\__     UPDATE delta={delta}')
                if dt < 0.:
                    # print('WARN: Too small dt, Rd too large or Ceff decreased too fast')
                    break
            if not nr_converged:
                print('WARN: smallNR not converged after many iters.')
                if Ceff_old is not None:
                    print('WARN: -> trying to lower Ceff step..')
                    Ceff = (Ceff + Ceff_old) / 2.
                    dt, t0 = dt_old, t0_old
                    continue
            
            # # prevent math.exp overflow
            # try:
            #     exp_term1 = np.exp(dt * poles / REGION)
            # except FloatingPointError:
            #     exp_term1 = np.exp(np.clip(dt * poles / REGION, -700, 700))
            # except OverflowError:
            #     exp_term1 = np.exp(np.clip(dt * poles / REGION, -700, 700))
            # try:
            #     exp_term2 = math.exp(-dt / (Ceff * Rd * REGION))
            # except OverflowError:
            #     exp_term2 = 1e-300
            # try:
            #     exp_term3 = math.exp(-dt / (Ceff * Rd * REGION))
            # except OverflowError:
            #     exp_term3 = 1e-300
            # Qpidt = np.dot(residues_mat[0], (-REGION + exp_term1 * REGION - dt * poles) / (dt * poles**2 * Rd * REGION))
            # QCeffdt = Ceff / REGION + Rd * Ceff**2 * (-1. + exp_term2) / dt
            # QCeff_quad_coeff = Rd * (-1. + exp_term3) / dt
            # inside_sqrt = 1. / REGION**2 + 4. * QCeff_quad_coeff * Qpidt
            # if inside_sqrt < 0:
            #     inside_sqrt = 0.0 
            # Ceff_new = (-1. / REGION + math.sqrt(inside_sqrt)) / (2. * QCeff_quad_coeff)
            
            
            Qpidt = np.dot(residues_mat[0], (-REGION + np.exp(dt * poles / REGION) * REGION - dt * poles) / (dt * poles**2 * Rd * REGION))
            QCeffdt = Ceff / REGION + Rd * Ceff**2 * (-1. + math.exp(-dt / (Ceff * Rd * REGION))) / dt
            QCeff_quad_coeff = Rd * (-1. + math.exp(-dt / (Ceff * Rd * REGION))) / dt
            Ceff_new = (-1. / REGION + math.sqrt(max(0., 1. / REGION**2 + 4. * QCeff_quad_coeff * Qpidt))) / (2. * QCeff_quad_coeff)
            if Ceff_new > Ceff:
                Ceff_new = (Ceff + Ceff_new) / 2.
            # print(f'| NEW Ceff={Ceff_new} (old: {Ceff})')
            # assert Ceff_new < Ceff, 'Should most likely be decreasing.'
            Ceff_old = Ceff
            Ceff = Ceff_new
            delay_old = delay
            slew_old = slew
            dt_old, t0_old = dt, t0
        else:
            print('WARN: Ceff and delay not converged after many iters.')

        return Ceff, delay, slew, dt, t0

    def calc_single_timepoint_nr(self, dt, node_id, v, init_t, residues_mat, poles):
        t = init_t
        last_t = t
        last_absf = 100.
        for i in range(20):
            f = self.calc_waveform(t, dt, residues_mat, poles, node_id=node_id) - v
            absf = math.fabs(f)
            g, _ = self.calc_waveform_grad(t, dt, residues_mat, poles, node_id=node_id)
            # print(f'single timepoint nr: i={i}, t={t}, f={f}, g={g}')
            if absf < 1e-3:
                # print(f'timepoint of node {node_id} v {v}: converged in {i} iters')
                break
            if not (absf < last_absf):
                # print('lowering step to solve the timepoint')
                t = (t + last_t * 3.) / 4.
                continue
            last_absf = absf
            new_t = t - f / g
            last_t = t
            t = new_t
        else:
            print('WARN: single timepoint not converged after many iters.')
        return t
    
    def Calc_Elmore(self, net, inslew, node):
        rc = self.Parasitics[net]
        if net not in self.Elmore.keys():
            assert len(rc.coupling_caps) == 0, \
                'Elmore delay don\'t support coupling caps'
            assert len(rc.grounded_caps) <= rc.n and len(rc.ress) == rc.n - 1, \
                'Must be a tree'
            caps = [0. for i in range(rc.n)]
            load = [0. for i in range(rc.n)]
            delay = [0. for i in range(rc.n)]
            ldelay = [0. for i in range(rc.n)]
            beta = [0. for i in range(rc.n)]
            sibling = [[] for i in range(rc.n)]
            for a, c in rc.grounded_caps:
                caps[a] = c
            for a, b, r in rc.ress:
                sibling[a].append((b, r))
                sibling[b].append((a, r))
            
            def dfs_load(u, f):
                load[u] = caps[u]
                for v, r in sibling[u]:
                    if v == f: continue
                    dfs_load(v, u)
                    load[u] += load[v]
                    
            def dfs_delay(u, f):
                for v, r in sibling[u]:
                    if v == f: continue
                    delay[v] = delay[u] + r * load[v]
                    dfs_delay(v, u)

            def dfs_ldelay(u, f):
                ldelay[u] = caps[u] * delay[u]
                for v, r in sibling[u]:
                    if v == f: continue
                    dfs_ldelay(v, u)
                    ldelay[u] += ldelay[v]
                    
            def dfs_beta(u, f):
                for v, r in sibling[u]:
                    if v == f: continue
                    beta[v] = beta[u] + r * ldelay[v]
                    dfs_beta(v, u)

            dfs_load(0, None)
            dfs_delay(0, None)
            dfs_ldelay(0, None)
            dfs_beta(0, None)
            
            self.Elmore[net] = (delay, beta)
        else:
            (delay, beta) = self.Elmore[net]
        
            
        outslew = []
        for d, b in zip(delay, beta):
            if inslew ** 2 + d ** 2 - 2 * b > 0:
                outslew.append(math.sqrt(inslew ** 2 + d ** 2 - 2 * b))
            else:
                outslew.append(inslew)
        
        return delay[rc.name2id[node]], outslew[rc.name2id[node]]
        
if __name__ == "__main__":
    Cal = Delay_Calculator()
    import TimingLib
    TB = TimingLib.TimingLib()
    Delay, Trans = TB.NLDM[('ssgnp', '0p63', '125')]['CKND1BWP16P90ULVT'].delay[('I->ZN', 'r', 'negative_unate','combinational', '*', '*')], \
        TB.NLDM[('ssgnp', '0p63', '125')]['CKND1BWP16P90ULVT'].trans[('I->ZN', 'r', 'negative_unate','combinational', '*', '*')]
    Cal.Calcu_Pin_Ceff_Delay_Trans('n20', 0.1, Delay, Trans)
    delay, beta = Cal.Calc_Elmore('n1410', 0.1, 'U754/A1')
    print(delay, beta)
    #Delay, Trans = TB.NLDM[('ssgnp', '0p63', '125')]['INR2D1BWP16P90ULVT'].delay[('CP->Q', 'r', 'non_unate','rising_edge', '*', '*')], \
    #    TB.NLDM[('ssgnp', '0p63', '125')]['INR2D1BWP16P90ULVT'].trans[('CP->Q', 'r', 'non_unate','rising_edge', '*', '*')]    